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Abstract. In the last decade, the perfectly matched layer (PML) approach has proved 
a flexible and accurate method for the simulation of waves in unbounded media. 
Most PML formulations, however, usually require wave equations stated in their stan- 
dard second-order form to be reformulated as first-order systems, thereby introducing 
many additional unknowns. To circumvent this cumbersome and somewhat expen- 
sive step, we instead propose a simple PML formulation directly for the wave equation 
in its second-order form. Inside the absorbing layer, our formulation requires only two 
auxiliary variables in two space dimensions and four auxiliary variables in three space 
dimensions; hence it is cheap to implement. Since our formulation requires no higher 
derivatives, it is also easily coupled with standard finite difference or finite element 
methods. Strong stability is proved while numerical examples in two and three space 
dimensions illustrate the accuracy and long time stability of our PML formulation. 
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1 Introduction 

The accurate and reliable simulation of wave propagations in unbounded media is of fun- 
damental importance in a wide range of applications. The perfectly matched layer (PML) 
approach [6J has proved a flexible and accurate method for the simulation of waves. It 
consists in surrounding the computational domain by an absorbing layer, which gener- 
ates no reflections at its interface with the computational domain; hence, it is perfectly 
matched. Inside the absorbing layer a damping term is added to the wave equation, 
which acts only in the direction perpendicular to the layer. This approach is analogous to 
the physical treatment of the walls of an anechoic chamber and provides an alternative 
to absorbing or nonrelfecting boundary conditions Blllll2l[T4l - [T6| . 

The initial PML formulation of Berenger |6| was based on splitting the electromagnetic 
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fields into two parts, the first containing the tangential derivatives and the second con- 
taining the normal derivatives; damping was then enforced only upon the normal com- 
ponent. Later Abarbanel and Gottlieb [IJ showed that Berenger's approach was only 
weakly well-posed due to the unphysical splitting of the field variables. Several strongly 
well-posed approaches have been suggested since, some of which were shown to be lin- 
early equivalent [|2l l20| . 

The PML approach has proved very successful in practice, because of its simplicity, ver- 
satility, and robust treatment of corners. Once discretized and truncated at a finite thick- 
ness, the layer is no longer perfectly absorbing and the optimal damping parameters need 
to be determined via numerical experiments. Stability properties of the PML approach 
has been analyzed in several works, such as in among others. 

The best implementation in the time domain is still under debate. Most PML formula- 
tions require wave equations stated in their standard second-order form to be reformu- 
lated as first-order hyperbolic systems, thereby introducing many additional unknowns. 
Here we propose instead a simple PML formulation directly for the second-order wave 
equation both in two and in three space dimensions. Our formulation also requires 
fewer auxiliary variables than previous formulations for the second-order wave equa- 
tion - see ll3ll5l[T9ll, for instance. 

Our paper is organized as follows. In Section 2 we derive a PML formulation for the 
wave equation in its standard second-order form. By judiciously choosing the auxiliary 
variables in the Laplace transformed domain, the resulting PML modified equations re- 
quire only two auxiliary variables in two dimensions and four auxiliary variables in three 
dimensions inside the absorbing layer. Next, in Section 3 we prove stability of our PML 
formulation by using standard theory from IITSll . The finite difference discretization of the 
PML modified wave equation is shown in Section 4. In Section 5, our numerical results 
both in two and three space dimensions demonstrate the accuracy and long time stability 
of the PML formulation. 

2 PML formulation 

We consider a time dependent wave field u propagating through unbounded three di- 
mensional space and assume that all sources and initial disturbances are confined to the 
rectangular domain 0= [— fli,fli] x [—a2,a2\ x [— fl3,fl3], ai,a2,a3>0. Outside O, we further 
assume the speed of propagation c> to be constant; hence, all waves are purely outgo- 
ing in the unbounded exterior R"'\n. Inside Q, the wave field u{xi,X2,X2,,t) satisfies 




(2.1) 
(2.2) 
(2.3) 



U = Uq 



Ut = Vo 



t = 0. 



We wish to truncate the unbounded exterior and thereby restrict the computation to the 
finite computational domain Q. In doing so, we need to ensure that all waves propagat- 
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ing outward leave Q without spurious reflection. Thus we shaU surround O by a per- 
fectly matched layer (PML) of thickness Li, i =1,2,3, in each coordinate which is designed 
to absorb the waves exiting Q. Inside the absorbing layer, u then satisfies a modified 
wave equation whose solutions decay exponentially fast with distance from the compu- 
tational domain. 

Following llilll, we let u denote the Laplace transform of u, defined as 

/■CO 

u = u(x,s)= / e^^u(x,t)dt, see. (2.4) 
Jo 

Outside Q, u then satisfies the Helmholtz equation. 

Next, we introduce the coordinate transformation 

Xii-^Xi:=Xi + - / ^i{x)dx, i = 1,2,3, (2.6) 



SJO 

where the damping profile ^, is positive inside the absorbing layer, | > Uj, i = 1,2,3, but 
vanishes inside O. If we now require u to satisfy the modified Helmholtz equation in 
those stretched coordinates, 

dxi V 9^1 J dx2 V dx2 J 9x3 V 9x3 y ' 

it is well-known that u will remain unaltered inside O, but decay exponentially fast inside 
the layer; hence the absorbing layer will be perfectly matched. In fact, the (unsplit) PML 
modified Helmholtz equation (|2.7|) in the Laplace transformed domain is standard [HQ. 
The difficulty lies in transforming (|2.7|l back to the time domain, without introducing 
high order derivatives or too many auxiliary variables. 

From (I2.6p , we observe that partial differentiation with respect to x, is related to partial 
differentiation with respect to the physical coordinate, x„ through 

3 s 3 

9x,- s + ^,3x,' ^^'^^ 
We now let ji = s), i = 1,2,3 denote 

7.:=1 + -- (2.9) 
s 

Then, by replacing partial derivatives according to (|2.8|l and multiplying the resulting 
expression by 717273, we rewrite (|2.7|) in physical coordinates as 

s 717273" dxi\ 7i 9xiy 9x2 V 72 9x2/ 9x3 V 73 9x3/ 
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From (|2.9|) we derive after some algebra the following identities: 

S7i (s+^i)s 



7371 



_l_ (^3+^1-^2)5+^3^1 ^2 11) 



572 (5+^2)5 
7172 ^-^ I (^1 + ^2-^3)5+^1^2 

573 (5+^3)5 

By using (|2ll]l in (|2l0)) we find 

s2 + s(Cl+^2 + C3) + (^1^2 + ^2C3 + C3^l) + ^^)" 

3xi V 9^1 y 9^2 V 9x2 / 9^3 V 9x3 J 
I 9 / 2/ (^2 + ^3-^i)5 + ^2^3 \ 9m\ 9 / . (^3+^1-^2)5+^3^1 \9i2 
9x1 V V (s+Ci)5 ;9xiy 9x2 V V (5 + ^2)5 /9x2 

I 9 / 2/ (gl+g2-^3)5 + glg2 \ du\ 
^9X3V V (s + ^3)5 ^9X3/ 

Next, we introduce the auxiliary functions ip and (p = {(pi,(p2,(p3)^ , 



(2.12) 





1. 

-u, 
s 






4>2 = 




h = 





or equivalently 



is + ^2)<p2 = c'{ (^3 + ^1-^2) 4 



^l)+ 


^2^3^ 
5 . 


\ du 


^2) + 


^3^1^ 
5 / 


\ du 


^3) + 


^1^2^ 
5 . 


\ 9w 
/9^' 
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Finally, we use the above relations in (|2.12|) and transform the resulting equations back to 
the time domain, which yields the PML modified wave equation 



Utt + {Cl+^2 + ^3)Ut + {^1^2 + ^2^3 + ^3^l)u = V-{c^Vu)+V-(p-Cl^2^3^p, 

(p^ = Ti(p + c^T2Vu+c^T3Vxp, 



(2.13) 



tpt = u, 



where 



-^1 





■ 




"^2 + ^3-^1 











-C2 





, r2 = 





^3 + ^1-^2 


















^1- 


Hi 





















^3^1 



















^1^2. 











and T3 ■ 



In the interior of O, the damping profiles ^,,i = 1,2,3 and the auxiliary variables ^, xp van- 
ish; hence, (|2.13|) reduces to (|2.1|) in Q. Because our PML formulation (|2.13|) requires only 
four auxiliary scalar variables (pi,(p2,(p3/^ inside the layer and no high order derivatives, 
its implementation is not only straightforward but also cheap to implement. 
In two space dimensions, ^3 and (p3 and xp vanish and our PML formulation reduces to 



Utt + {^i + ^2)ut + ^H2U = V ■ {c^Vu)+V -cp, 



(|>^ = Tllp+c T2VU, 



(2.14) 



where 



Hi 






Hi 



r2 



Cl-C2 



Remarkably only two auxiliary functions are needed here. 

The choice of the damping profiles ^/(x) > 0, i = 1,2,3 is arbitrary; it can be constant, 
linear, or quadratic among others. In our computations, we always use 




\xj—aj\ 



2n 



for \xi \ <ai, i = 1,2,3 
for fl,- < I X, I < Uj + Lj, i - 



: 1,2,3. 



(2.15) 



Because ^i{x) is twice continuously differentiable throughout the interface at |x;| = fl,, 
no special transmission conditions are needed there. The constant ^, depends on the 
discretization and the thickness of the layer, which in practice is truncated by a homo- 
geneous Dirichlet (or Neumann) boundary condition. Then the relative reflection, R, is 
given by 

^, = |:log(^), / = 1,2,3. (2.16) 
In Figure 1 we show damping profiles for different values of ^j. 
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Figure 1; The damping profile ^/(x,) given by (|2.15|l is shown for different values of ^j-, with c = l and L; = 0.1 . 

3 Stability 

We now establish the stability and well-posedness of our PML formulation, first in two 
and then in three space dimensions, where we assume that the absorbing layer extends 
to infinity. Here we follow standard stability theory for hyperbolic systems ||T8| . which 
we briefly recall below. 
Consider a general Cauchy problem, 

Ut = p(^^^U, 0<t<T, UeW, (3.1) 

where P{dx) denotes a linear differential operator, with initial conditions 

U{x,0) = Uo{x), xeM.^. (3.2) 

Following ItTsU , the Cauchy Problem is weakly (resp. strongly) well-posed, if the solution 
U{-,t) satisfies 

\\U{-,t)\\L,<Ke^'\\U{-,0)\\H^ (3.3) 

with s > (resp. s = 0). The Cauchy Problem is weakly (resp. strongly) stable, if the solution 
U{-,t) satisfies 

\\U{-,t)\\L,<K{l + ty\\U{-,0)\\H^ (3.4) 

with s > (resp. s = 0). A necessary and sufficient condition for weak well-posedness (resp. 
stability) is that all eigenvalues A of the operator P{ik)) satisfy 

K{A(P(*))}<C, fcGR, (3.5) 

with C > (resp. C = 0) independent of k. For strong well-posedness (resp. stability), the 
corresponding eigenvectors must also be complete. 

By rewriting the PML-modified wave equations (|2.13p , (|2.14|) as a first-order hyperbolic 
system and applying the stability theory from IITSlI delineated above, we can prove the 
following two stability results. 
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Theorem 3.1. The Cauchy problem for the PML formulation (|2.14|) in two space dimen- 
sions is strongly stable for ^1,^2 > 0. 
proof) 

For simplicity, we assume that ^1,^2 are constant; note, however, that the stability theory 
from IjlSf extends to smoothly varying coefficients. We introduce the new variable v to 
rewrite the first equation in (|2.14|) equivalently as 



Mt = — ^2W + divv, vt = — ^iv+c Vm + ^. 
By using (I3.6D . we now rewrite (|2.14P as a first order hyperbolic system: 

Ut = AUx + BUy + C, 

where 



A: 



and C 





Ut = 


{u,(p 


1'^ 


f'2,Vi,V2) 


1 

r 
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0' 
















1 




) 




































, B = 

















C2 


























































X2 





0" 
















^1 

























^2 






















^1 








































(3.6) 

(3.7) 
(3.8) 



(3.9) 



By using a symbolic algebra program we find that the eigenvalues of the principal part 
of P(zfc) for (|321lare 

X{V{ik)) = ±ic{k\^'k\f'^. (3.10) 



Thus, 



3?{A(P(zfc))} = 0, 



(3.11) 



while the corresponding eigenvectors are also complete for all ^1,^2 > 0. Therefore, since 
C is a diagonal matrix with negative entries for ^1,^2 > 0, we conclude that (|2.14|) is 
strongly stable. 



Theorem 3.2. The Cauchy problem for the PML formulation (|2.13|l in three space dimen- 
sions is strongly stable, if at least two = 0,; = 1,2,3, and weakly stable, otherwise, 
proof) 

We introduce the new variable v to rewrite the first equation in (12.1311 as 

Wf = -^2W + divv-^3t/?, vt = -^iv+c2VM + </>. (3.12) 
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By using (|3.12|l , we can rewrite (|2.13|) as a first order hyperbolic system: 



where 



A: 



and 



C- 



Ut = AU^ + BUy + CU^ + D, 



Ut = {u,(pi, (p2, (p3, Vi, V2, V3, 1/?) ^, 







1 



c2(^2 + ^3-^i) 


























. 

1 " 



^3^1 









. 

1 " 





c2(^l+^2-^3) ClC2 



c2(^3 + ^l-^2) 



c2 





















(3.13) 



(3.14) 



(3.15) 



(3.16) 



(3.17) 



By using a symbolic algebra program, we find that the eigenvalues A of P{ik) for (|3.13|) 
are 

A{P{ik)) = ±ic{kl+kl+kl)^/^. (3.18) 



Thus, 



^{A{Piik))} = 0, 



(3.19) 



while the corresponding eigenvectors are also complete, if at least two Q = 0,] = 1,2,3; 
else, they are not complete. Therefore, since D is a diagonal matrix with negative entries 
for ^1,^2,^3 > 0, we conclude that (|2.13l) is strongly stable, if at least two Q=0, and weakly 
stable, otherwise. 



4 Finite difference discretization 



Here we show how to discretize (|2.13|) with standard second-order finite differences on a 
uniform mesh at grid points ,-^=Xi o+z^Ax/, with i=l,2,3 and i(=0,l,...,M(. For the time 
discretization we use a constant step size At and denote the time levels by t„ = to + nAt, 
n=0,l,...,N. Inside the absorbing layer, we further introduce a space-time staggered grid 
at locations x--^^i = Xifi + [ii + j) Axj, i = 1,2,3 and times t^^i = fo + (n + 2) Af. Then the 

numerical solution u^^j^, which approximates u at grid point (:^i,/,^2,;,^3,*:) and time 
satisfies 



— — + {^li + ^2j + ^3k) + {^li^2j + ^2j^3k + ^3k^li) ulj,k 



i+lj,k^'i+l,j,k (^i+lj,k + ^i-l,j,k'^ ^^li,k + '^i- lj,k^'i-l,j,k 

A^ 

^ '^Ij+lkKj+U - i^lj+lk^^lj-lk^Kj,k + '^lj-lk^li-U 



AX2^ 

_^ ^lj,k+ 1 Ki,k+1 ~ ^'^lj,k+ \ + '^IjM- \ ^ "^^7'^ ^ '^lj,k- \ Kj,k-l 



AX3^ 

Xn Xn Xn Xn Xn Xn ,"+7 , 

Jli+Uk'^^li-l^M ,'^2^+\,k-^2i4-\,k ^^3^^^^^^ ^ ^ ^ t,j,k +Hk 

+ A^^ + A^2 ^ A^^^ 2 



where the cell averages of the auxiliary functions (pi, (p2 and (ps are defined as 
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^2i,j+i,k = l {^li-\,j+\,k-\ +'Pli-ij+ik+\ +(p2i+y+ik-i +(Pli+i,j+i,k+\ 

^3i,j,k+\ = 4 {"^li-li-lk+l +'Pli-\,j+\,k+\ '^'^3i+\,i-\,k+\ ^'^3r+\,j+\,k+\, 



Concurrently with the above discretized wave equation, we also advance the (scalar) 
auxiliary variables t/^, ^j, j = 1,2,3 inside the absorbing layer by using standard finite 
differences. For xp we use 
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whereas for (pi we use 



li+^,i+l,k+l 



At 



^"u+lj+l,k+ 1 + 'f^li+ \ ,j+ \ M 



'2;- 



where 



i+^,j+^,k+^ 



^1 




1 , 1 — w' 






Axi 
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i,l+l,k^ 



u" 

i+l,i+\,k+\ 



i,j+\,k+\ 



Axi 



Here, the cell averages of u and \p are defined as 



i,j+\,k+\ 



^r,i+\,k+\ 



U,j,k+1 "T ^i,i+l,k "T "/,;+l,7c+l 



,.n+2 



11+2 



t,j,k +t,j,k+l + t,i+lk + t,j 



The finite difference approximations for ^2 and (|)3 are analogous. 



5 Numerical experiments 

Here we present numerical experiments that illustrate the accuracy, versatility and long- 
time stability of our PML formulation discretized with standard finite differences as in 
Section 4. In all cases we choose = 80 in the damping profile, which yields a relative 
reflection R ^ 10^^ for the the typical values c = 1 and L, = 0.1. At the exterior boundary 
of the absorbing layer we impose homogeneous Dirichlet boundary conditions. 



5.1 Point source in 2D 

First, we consider the wave equation (|2.H) in two space dimensions with constant speed 
of propagation c=l and zero initial conditions, Uq=Vq=0. The source term / corresponds 
to a truncated first derivative of a Gaussian: 



f{x,y,t)=S{x)3{y)h{t) 



(5.1) 
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with 

h{t) = j-^{e-^'(f°^-'^'), /o = 10Hz. (5.2) 
The grid spacing is uniform in Xi and X2, with Ax = 0.002. 

In Figure 2 we display snapshots of the numerical solutions at different times in Q = 
[—0.5,0.5]^, surrounded by a PML of width L = 0.1. We observe how the circular wave 
propagates outward essentially without spurious reflection from the PML. By time t = 
1 the wave has essentially left the computational domain. To assess the error in the 
numerical solution, we compute a reference solution in a much larger domain of size 
[—5.5,5.5] X [5.5,5.5], so that boundary effects are postponed to later times inside O. In 
Figure 3, the time evolution of the L^-error is shown for different values of the damping 
coefficient ^j. Until f = 8 we observe a steady decrease of the error over seven orders of 
magnitude, regardless of the value of ^„ which demonstrates the long-time stability of 
our method. Moreover, our formulation appears robust with respect to the parameter 
value [i. 



1 ^ 0.40729 ,10' 1=0,58548 .,o' t = 0.64912 




Figure 2: Point source in 2D: snapshots of the numerical solutions at different times in Q = [—0.5,0.5]^, 
surrounded by a PML of width L = 0.1. 



5.2 Heterogeneous medium in 2D 

Next, to illustrate the versatility of our PML formulation, we consider the homogeneous 
wave equation l|2.ip in a heterogeneous medium with varying wave speed c=c(x2), given 
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Figure 3: Point source in 2D: time evolution of the L2-error for different damping coefficients 



by 



c{Xi,X2)=< 



0.5, ifx2<-b 
l + l^ + 2^sin(^), if|x2|<fo 
1.5, otherwise. 



(5.3) 



We set b = 0.95 which yields the vertical velocity profile shown in Figure IH The initial 
conditions are 

u\t=o = Uo{xi,X2) and Ut\t=o = 0, (5.4) 



where 

Mo(Xi,X2) 



(4(xi +0.4) (0.4-:«i) ) sin(37rx2), if -0.4 <xi< 0.4,-1 <X2<1 



0, otherwise. 
Here O is the square domain [—1,1] x [—1,1], surrounded by a PML of width L=0.2. The 



1.5 
1.25 

CO 

->S 1 



0.75 



0.5 



-0.5 0.5 



Figure 4: Heterogeneous medium in 2D: varying wave speed c given by (|5.3 
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finite difference grid is uniform with grid spacing Ax = 0.004. In Fig|5l we display snap- 
shots of the solution at different times, where again the last frame is purposely chosen 
at a much later time. In spite of the varying wave speed and the glancing angle of in- 
cidence along the vertical artificial boundaries, the waves are damped without spurious 
reflection. Even at much later times we do not observe any instability in the numerical 
scheme. 



„j 1 = 0.10182 1 = 0.50912 




Figure 5: Heterogeous medium in 2D: snapshots of the numerical solution are shown at different times in 
n=[-l,l]2, surrounded by a PML of width L = 0.2. 



5.3 Point source in 3D 

Finally, we consider the wave equation (|2.1|) in three space dimensions with zero initial 
conditions and the same point source / as in (|5.1|l . The grid spacing is uniform in x-[, X2 
and Xj, with Ax = 0.006. In Figure 6, we display snapshots of the numerical solutions at 
different times in O = [—0.5,0.5]^, surrounded by a PML of width L = 0.1. We observe 
how the spherical wave propagates outward essentially without spurious reflection from 
the PML. By time t = l the wave has essentially left the computational domain. Again we 
observe no instabilities in the numerical solution even at much later times. 

6 Concluding remarks 

We have presented a PML formulation for the wave equation in its standard second-order 
form. It distinguishes itself from known formulations by its simplicity and the small 




Figure 6: Point source in 3D: snapshots of the numerical solution are shown at different times in 0= [—0.5,0.5] 
surrounded by a PML of width L = 0.1. 
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number of auxiliary variables needed inside the absorbing layer. We have proved that 
the continuous Cauchy problem with the unbounded PML is stable and well-posed. Our 
numerical results in two and in three space dimensions with standard finite differences 
illustrate the accuracy, versatility and long-time stability of our PML formulation. 
Because it involves no high space or time derivatives, our PML formulation easily fits 
continuous or discontinuous Galerkin formulation for use with finite element methods 
|T0l[T3ll . It also immediately generalizes to Maxwell's equations in second-order form. 
Current work involves the extension to second-order wave equations in complex elastic 
and poro-elastic media, and will be reported elsewhere in the near future. 
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